A Rigorously Justified Robust Algebraic Preconditioner for High-contrast Diffusion Problems

In this paper we analyse the robustness of (algebraic) multigrid preconditioners applied to linear systems arising from finite element approximations of elliptic PDEs with high-contrast coefficients. Problems with high-contrast coefficients are ubiquitous in porous media flow applications. Consequently, development of efficient solvers for high-contrast heterogeneous media has been an active area of research. Here we are particularly concerned with the convergence of a family of algebraic preconditioners that exploit the binary character of high-contrast coefficients (see also [1]). We consider preconditioners for piecewise linear finite element discretisations of boundary-value problems for the model elliptic problem −∇ · (α∇u) = f , (1) in a bounded polygonal or polyhedral domain Ω ⊂ Rd, d = 2 or 3 with suitable boundary conditions on the boundary ∂Ω. The coefficient α(x) may vary over many orders of magnitude in an unstructured way on Ω. A finite element discretisation of (1) yields the linear system Au = f , and it is well-known that the conditioning of A worsens when the mesh T h is refined or when the heterogeneity (characterised by the range of α) becomes large. It is of interest to find solvers for this system which are robust to the heterogeneity as well as to the mesh width h. In the literature there are many papers devoted to the efficient solution of this problem with a rigorous justification when discontinuities in α are simple interfaces which can be resolved by a coarse mesh. Even if suitable coefficient-resolving coarse meshes are not available, good performance of Krylov-based methods can still be achieved by standard preconditioners when there is a small number of unresolved interfaces. This is because the preconditioning produces a highly clustered spectrum with correspondingly few near-zero eigenvalues ([2,5]). For more general complicated heterogeneous high-contrast media, recent progress was made in [3,4]. At the same time it is well-known that algebraic multigrid (AMG) procedures also produce optimal robust solvers for such heterogeneous problems, but so far theoretical justification of their robustness with respect to coefficient variation is lacking. The family of algebraic preconditioners proposed in this paper can be constructed in similar ways as AMG preconditioners by identifying strong and weak couplings in the stiffness matrix. However for this new family we can prove the robustness and we demonstrate this on a sequence of model problems. Moreover, our numerical experiments show that for sufficiently high contrast the performance of our new preconditioners is almost identical to that of the Ruge and Stüben AMG preconditioner, both in terms of iteration count and CPU-time ([1]). To give some more details, the first (algebraic) phase of our family of preconditioners involves partitioning of the degrees of freedom into a set corresponding to a “high-permeability” region ΩH and a “low-permeability” region ΩL. Note that (for sufficiently high contrast) this can easily be obtained by using a strong-connection criterion similar to that used in AMG algorithms. Thus the stiffness matrix can be partitioned into A = [ AHH AHL ALH ALL ] . After a little algebra, the exact inverse of A can be written: A−1 = [ I −A−1 HHAHL 0 I ][ A−1 HH 0 0 S−1 ][ I 0 −ALHA HH I ] (2) where S = ALL −ALHA HHAHL is the Schur complement of AHH in A. A singular perturbation analysis can now be devised to explain the properties of the subblocks in (2) and to derive sharp bounds for the eigenvalues of the preconditioned stiffness matrix. Arguments of this type were first used in the context of condition number analyses for additive Schwarz methods in [2]. More recently this approach was refined to treat the more complicated problem of analysing multigrid preconditioners in [6]. Here we use the singular perturbation-type analysis in a different context. Suppose for simplicity that α|ΩH = α̂ 1 and that α|ΩL = 1. (Note however, that our method and our analysis are not restricted to this piecewise constant model situation.) It is clear that α̂AHH = NHH +O(α̂−1) , as α̂→∞, where NHH is the matrix corresponding to the pure Neumann problem for the Laplace operator on ΩH . This shows that (after scaling by α̂−1) AHH can be preconditioned robustly and efficiently by standard multilevel methods, such as geometric multigrid, with a performance independent of h and α̂. Moreover the analysis of AHH as α̂→∞ has important implications for the behaviour of S. We show that in this case S = S(∞) + O(α̂−1) , where S(∞) is a low rank perturbation of ALL. The rank of the perturbation depends on the number of disconnected components (“islands”) in ΩH . This special limiting form of S allows us to build robust approximations of S−1, e.g. combining solves with ALL (again available robustly using standard multilevel methods) with the Sherman-Morrison-Woodbury formula. Finally, we show that the application of the two remaining blocks in (2) corresponds (in the high-contrast limit) to deflation with respect to certain low frequency eigenvectors.